DAHLGREN DIVISION 
NAVAL SURFACE WARFARE CENTER 


Dahlgren, Virginia 22448-5100 


NSWCDD/TR-05/89 


AN INVERSE OF THE CIRCULAR COVERAGE 
FUNCTION 


BY ARMIDO DIDONATO 


FORCE WARFARE SYSTEMS DEPARTMENT 


JANUARY 2005 


Approved for public release; distribution is unlimited. 


REPORT DOCUMENTATION PAGE 


Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, search existing data sources, gathering and 
maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden or any other aspect of this collection of information, including 
suggestions for reducing this burden, to Washington Headquarters Services. Directorate for information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302, 
and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188), Washington, DC 20503. 


1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED 
January 2005 Final 


4. TITLE AND SUBTITLE 5. FUNDING NUMBERS 


An Inverse of the Circular Coverage Function 


6. AUTHOR(s) 
Armido DiDonato 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION 


Commander REPORT NUMBER 


Naval Surface Warfare Center NSWCDD/TR-05/89 
Dahlgren Division (Code T10) 

17320 Dahlgren Road 

Dahlgren, VA 22448-5100 


9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 


11. SUPPLEMENTARY NOTES 


12a. DISTRIBUTION/AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE 


Approved for public release; distribution is unlimited. 


13. ABSTRACT (Maximum 200 words) 


This report describes an algorithm, INCIR, to evaluate an inverse of the Circular Coverage Function, CIR. Given a 
two-dimensional normal distribution function with mean zero and equal unit standard deviations and a target located at 
(h, k); R, the radius of the circle, C, with center at (h, k) is determined for a specified probability of a shot falling 
within C. 


A Fortran 77 subroutine, INVCIR, is available that is based on INCIR, which produces R to eight significant digits for 
10°<P<1-10%and0<d=(h*+k’)'’?< 10°. 


A table of R versus P and d is included. 


14. SUBJECT TERMS 15. NUMBER OF PAGES 
algorithm, INCIR, inverse, Circular Coverage Function, CIR, INVCIR 48 


16. PRICE CODE 


17. SECURITY CLASSIFICATION | 18. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION OF 


OF REPORTS OF THIS PAGE OF ABSTRACT ABSTRACT 


UNCLASSIFIED UNCLASSIFIED UNCLASSIFIED UL 


tandard Form 298 (Rev 2-8 
Prescribed by ANSI std. Z39-18 
298-102 


NSWCDD/TR-05/89 


BLANK PAGE 


NSWCDD/TR-05/89 


FOREWORD 


The algorithm described in this report is the basis for the Fortran software of an 
important statistical function that is not contained in the NSWC Library of Mathematics 
Subroutines. The software can be used in targeting studies when statistical confidence 
regions are required. 

Dr. John Crigler (B10) supplied the external distribution list. 

The editorial assistance of David Bozicevich (B60) is appreciated. 

This document was reviewed by Robert G. Hill, Head, Warfare Systems Division. 


Approved by: 


Rul. \ then 


REUBEN S. PITTS, Head 
Force Warfare Systems Department 


ill 


NSWCDD/TR—05/89 


BLANK PAGE 


Section 
I 
II 
III 
IV 


NSWCDD/TR—05/89 


CONTENTS 


Page 
INTRODUCTION sea isin cit ag eter esate a bon nav aduate aah nay ae ah 1 
ANALYSIS TO DETERMINE R (R= Ru)..........0.0c0ceceeeeceseeeeenes 2 
DESCRIPTION OF THE FORTRAN 77 SUBROUTINE INVCIR........... 4 
RE EER ENC BS 2 coehics soi tei wie pte ciarabarnann Yaoi aera saad anes don Reyes 5 
APPENDIX A—TABULATION OF R AS A FUNCTION OF P and d...... A-l 
DISTRIBU BION. issn hohe 0h ates, utnie « geicanha Na Be hla ecussiiins (fad (1) 


NSWCDD/TR—05/89 


BLANK PAGE 


vi 


NSWCDD/TR—05/89 


I. INTRODUCTION 


The Circular Coverage Function, CIR, defines a probability function P(R,d), [2]. CIR 
gives the probability that a shot falling under a bivariate uncorrelated normal distribution 
with mean zero and equal standard deviations, u = v, lands in a circle, in the x—y plane, of 
radius R centered at x = h, y =k. It can be expressed as 


R?—(x—h)? A oe 


h+R 
POR, B/a, Bf) = 55 exp {slr + (Pp dyads 1) 


Transforming to polar coordinates, with x —h = rcos@, y—k=rsin@ yields 


Be. exp{—sl (eats (ATES OY rac (2) 


u 


P(R/u, h/u, k/u) = 


Qru? 


Since u = v, the distribution is circular, hence the circle of integration can be taken as offset 
along the positive x-axis. Thus (2) reduces to 


R 
P(R, d) = exp(—d?/2) / exp[—r?/2) Ip(rd) r dr, (3) 
where is 
Ip(w) = - | exp(w cos 6) dé. (4) 


Here I,(w) denotes the modified Bessel function of the first kind and zero order [1, p. 375] 
and 
R=R/u, h=h/u, k=k/u, d=+/(h? +k). (5) 


Introducing the MacLaurin series and the asymptotic expansion for Io, [1, p. 375, p. 377], 
recurrence relations are developed which are used to generate an algorithm that efficiently 
and accurately evaluates P (see [2}). 

The objective of this report is to describe an algorithm upon which a Fortran computer 
program, INVCIR, is based to determine R, given P,d. This program can be used as a tool 
in evaluating the accuracy of a weapon. 

For example, consider a target located at the origin of the x—y plane and a set of miss 
distances normally distributed in the x and y directions [6], [7] that are independent, with 
mean zero and standard deviations u, u. Then, approximating the standard deviation from 
the data, say ti, one asks what is the radius, R, of the circle, C, centered at (h,k), which 
contains P of the distribution, i.e., for which a shot has a probability P of falling within C 
[6], [7]. 

In Appendix A, the six-decimal-digit inverse P(R, d) table of [2] is extended from six-seven 
significant digits to seven-eight significant digits. 
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II. ANALYSIS TO DETERMINE R (R = Ru) 


In INVCIR, a special case is treated first. When d = 0, then P = 1 — exp(—R?/2). 
Therefore, 


R = \/—2log (1 — P), if d=0, (IR = 0)!, (6) 


In order to find R in the general case, where d # 0 and P = P are given, the classi- 
cal Newton-Raphson root-finding procedure (N—R) is used [4, p. 129]. Denoting the nth 
approximation or iterate for R by Ry, with d constant, 


exp(—d?/2) ie exp[—r?/2) Ip(rd) rdr — P 


Ratt = Ra - (OP /OR) nex, 6. t= One est) (7) 
where a 24 
-- 
ae = Rexw(-F“)h(Ra (8) 
or if R > 1.7 and Rd > 16 (see [2]) then 
OP | (R —d)? 
m= Rexp(-S 57) fexp(—Ra) fo(R 4) (9) 


There are two basic problems associated with using N-R: 


1. A condition for stopping the iteration procedure is needed, i.e., deciding when an ac- 
ceptable approximation for R has been reached. 


2. A value for Ro, which initiates the N—R, is needed so that rapid convergence is assured. 


The first problem was resolved by requiring that the iteration be stopped when 


| Ra a a | es Rn41 eps, (10) 


where eps is a prespecified small positive number. It was found unnecessary to utilize the 
stopping rule: 
|P(Ray1,d) — P| < Peps. (11) 


The second problem required extensive testing and some analysis to resolve satisfactorily. 
The solution adopted distinguishes this work from that carried out in [2]. The reasoning will 
be heuristic. 

Values for Ro were considered using the results of [3], but those did not do as well, in 
some cases, as the starting values given below. 

By (6), it can be assumed d > 0. Then five choices were made for Rg depending on the 
input values of P and d. 


1TR is used in INCIR and defined in the next section. 
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First, consider expanding the integrand of (3) in a MacLaurin series in the integration 
variable r and integrating the result. Then one obtains, letting T = P exp (d?/2), 


Pay. a Ro (12) 


o. 
II 
mn 


a, = 1/2, ag = d?/16—1/8, ag = d*/384 — d?/48 + 1/48. 


Hence, Ro can be obtained, for a region of the P—d plane, by inverting the series. Thus, 
with y = V2T, D2 = d?/2, the first three terms of the inverted series are used. 


Ro = y [1+ (1 — D2) y?/8 + (17 D2? — 26 D2 + 13) y*/384], if T<.099. (IR=1) (13) 


Since, as noted above, 
P(R, 0) = 1 — exp(—R?/2), (14) 


we extend the use of (6) so that, if the inequality in (13) is not satisfied, then 
Ro = \/ —2log(1— P), if d<0.20. (IR = 2) (15) 


If the inequalities in (13) and (15) do not hold, the three following estimates are used for 
Ro with the region of the P—d plane where they apply specified by inequalities. If 


d > 2.5 or (P > .99 and d > .70), (16) 


then 
Ro = V(d+x)?+1 if d+x> .50, (iR=3) (17) 


where x is given by [1, p. 931, (26.2.2)] 
P= 1/vin f exp(—t?/2) dt, (18) 


and is obtained from the subroutine DPNI (P,Q, P—.5, x, IN)? of [5]. 
If the conditions in (16) hold, and d+x <.5, then [1, p. 940, (26.3.26)| 


Ro = (2 + a2) B3/1.50, 3 B =x/V9a +1 —1/(9a), (IR = 4) (19) 
where a is defined in (21). Finally, if the conditions in (16) do not hold, then from [3] 
Ro = y [22 (14+ d?)/(1+ D2)], (IR:= 5) (20) 
where z is given by [{1, p. 260, (6.5.1)] 
P= (a) 2 exp(—t) t*~+ dt, a = (1+ D2)*/(1+ d?), (21) 


2Q=1-P; IN is an error parameter. 
3The factor 1.50 was an adjustment of ours. 


NSWCDD/TR—05/89 


and is obtained from the subroutine GAMINV (a, z,zo, P,Q, IND) of [5]. zo is an optional 
initial approximation for z and IND reports the status of the results. 

IR is an output parameter of INVCIR that identifies which Ro was chosen. (See the next 
section.) 


III. DESCRIPTION OF THE FORTRAN 77 SUBROUTINE INVCIR 
In this section, the subroutine INVCIR is discussed further. Its call line is: 
Call INVCIR(P, d, R, ij, Ro, IR) 


This routine, with 28 supporting routines taken from [5], is written in Fortran 77 in double- 
precision. Thus, on an IBM PC approximately 15 decimal digits are available. The inputs 
are P and d. The outputs are R, ij, Ro, IR. IR plays a dual role, namely as an error parameter 
and also to indicate the choice for Ro as noted in the previous section. 
If IR is negative, then an unacceptable input has been given. R = —10!° is returned. In 
particular: 
If IR=—-—1 then P>lorP <0. 
If IR=-—2 then dO); 
The integer ij represents the number of calls to the N—R procedure. The output R = R/u. 
INVCIR is designed to give a minimum of eight significant digits for R when eps is set to 
5-10~-°, and P and d are constrained to 


10 <P<is 10%, 0<d< 108. 


The routine is also efficient. Extensive testing showed a maximum for ij of six and an average 
for ij of three, indicating that the choices for Ro are quite good. 

As an example of how to use the routine, assume P = .40, h = 3.0,k = 4.0,u =v = 
2.0. Then the input is P = .40,d = Vh?+k?/u = 2.50. INVCIR yields ij = 3, Ro 
2.45915620, R = 2.45457329. (see page A-9.) Hence R = R * u = 4.90914658. 
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APPENDIX A 


TABULATION OF R AS A FUNCTION OF P AND d 
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TABLE A-1. TABULATION OF R AS A FUNCTION OF P AND d 


Constant value of P on a row. Constant value of d on a column. 


Example 1: P= .31, d=.40, R=0.8963486 see Page A-5 
Example 2: P= .9990, d= 20.0, R= 23.113476 See Page A-26 
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